Viscoelastic Fracture of Biological Composites 



Eran Bouchbinder 1 and Efim A. Brener 1,2 

Chemical Physics Department, Weizmann Institute of Science, Rehovot 16100, Israel 
2 Peter Griinberg Institut, Forschungszentrum Jiilich, Jiilich 52425 Germany 



Abstract 

Soft constituent materials endow biological composites, such as bone, dentin and nacre, with viscoelastic 
properties that may play an important role in their remarkable fracture resistance. In this paper we calculate 
the scaling properties of the quasi-static energy release rate and the viscoelastic contribution to the fracture 
energy of various biological composites, using both perturbative and non-perturbative approaches. We con- 
sider coarse-grained descriptions of three types of anisotropic structures: (i) Liquid-crystal-like composites 
(ii) Stratified composites (iii) Staggered composites, for different crack orientations. In addition, we briefly 
discuss the implications of anisotropy for fracture criteria. Our analysis highlights the dominant lengthscales 
and scaling properties of viscoelastic fracture of biological composites. It may be useful for evaluating crack 
velocity toughening effects and structure-dissipation relations in these materials. 
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1. Introduction 



Biological composites, 


such as bone, dentin and nacre, exhibit outstanding mechanical properties ( 


Arzt et al. 


2003, 


Peterlik et al.. 


2006 


Fratzl and Weinkamer. 2007: Ritchie et al. 20091: Dunlop and Fratzl, 2010 


; Ritchie . 


201G 


Launev et al.. 


2010[ 


Ji and Gaol 2010). In Darticular. thev combine elastic stiffness and fracture resis- 



tance that is not yet achieved by synthetic composites of similar composition. Therefore, extensive recent 
efforts have been aimed at exploring the basic principles underlying the heterogeneous structures and defor- 
mation mechanisms of these materials, seeking guidelines for the development of novel synthetic composite s 
(|Fratzll . l2007t iMunch et all . 120081 lAntonietti and Fratzll . l2010i: IPunlop and Fratzll l2010i: I Ji and Gaol |2010| ). 



At the nano-scale, many biological composites consist of hard, plate-like, mineral crystals embedded in a 
soft protein matrix. The scale and spatial arrangement of the plate-like mineral crystals are believed to play 
a crucial role in endowing biological composites with their remarkable mechanical properties. For example, 
the nanometric dimensions of the mineral crystals in bone-l ike composites ha ve been proposed to be funda- 



mentally linked to the fracture resistance of these materials (jArzt et all 120031) . Furthermore, such biological 
nano-composites where shown to exhibit a generic nano-structure in wh ich the hard mineral crystals are 
arranged in a parallel staggered pattern ins ide the soft protein matrix (<Fratzj l2007t IPunlop and Fratzll . 
20101 Launev et al. . 2010t Ji and Gao . 201dh . Other nano-co mposites such as nacre, co mposed of parallel 
stratified arrays of hard mineral crystals, are of great interest (jPunlop and Fratzl . 2010l ). Our focus in this 
paper is on the macroscopic implications of these nano-structures and their constitutive behaviors. 

Recent effort has been devoted to the experimental characterization and modeling of the dependence of 
the fracture resistance of bi ological composites on the cra c k length, the so-call "crack extension resistance 
curve" (R-curve) behavio r (iNalla et~afl I2003L 12004 120051 iKinnev and Ritchie! . 120051: iKoester et all . 12008 : 
Launev and Ritchie . 20091 Launev et al. . 2010h . To the best of our knowledge, much less attention has been 
given to the increase in fracture resistance due to the finite velocity of cracks, in spite of the fact that 



even small velo c ities may generate a non-negligible contr i butio n to the fracture resistance (| Sasaki et al 
19931 lOkumural l2002al l2003l I Ji and Gaol l2004Hlvo eta!] . 120041) . In fact, the difference between the frac- 
ture toughness of hydrated and dehydrated biological composites may provide indirect evidence in favor 
of this possibility (jKruzid . I2003T ). This finite velocity effect can be attributed to the viscous component of 
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the mechanical respons e of the soft constituent materials in biological composites ([Puxkandl et all 12002 



Hazenberg et all [2006). Our goal in this paper is to explore this possible toughening mechanism in the 



framework of a single timescale linear viscoelastic model for various nano-structures. We also consider 
fracture initiation and some related anisotropic effects. 

To quantify fracture resistance, consider a crack moving at a velocity v and write the total energy release 
rate Gtot(v), i-e. the amount of energy needed to create a unit of crack surfaces, as 

Gtot{v) =G Q + G ms (v) , (1) 

where Go is the velocity-independent energy release rate and G V i S (v) is the finite-?; dissipation associated 
with viscous deformation. Note that for v > the energy release rate equals the fracture energy, T tot (v), 
which is a materials property, and hence Go should be identified with the critical energy release rate G c , i.e 
the fracture energy at the initiation of crack propagation. Nevertheless, we prefer to use the notation Go in 
Eq. (TTJ) since it allows interpreting our results even when Go<G c , i.e. under sub-critical conditions, where 
Go is the elastic energy release rate associated with a virtual incremental extension of the crack. Since 
G„i a (0)=0, we can expand G V i S (v) to the lowest order in v as 

G ms (v) ~ G ^w{£/d c ) + 0(v 2 ) , (2) 

where r is a typical viscous relaxation timescale, d c is the smallest scale cutoff for a continuum description 
in a given problem, £ is a quantity of the dimension of length and w(-) is a dimensionless function. Note that 
even though biological composites may exhibit a hierarchy of viscous relaxation times, throughout this paper 
we adopt the simplifying assumption that there exists only one relevant viscous relaxation time, associated 
with the soft constituent materials in biological composites. Generalization to several viscous relaxation 
times is rather straightforward. 

The generic two-dimensional fracture problem we consider in this paper consists of a crack of linear size 
L within a large body (i.e. a body whose linear dimensions are much larger than L) and under an applied 
remote tensile stress a^o that tends to open it. Dimensional analysis implies that we can write Go as 

Go = ^g(t/L) , (3) 

where E is a relevant elastic modulus, <?(•) is a dimensionless function and £ is a lengthscale. As implied 
above, the dimension of Go is energy per unit area. 

The ultimate goal of this paper is to calculate G V i s , either perturbatively as in Eq. Q or non- 
perturbatively as in Eq. (TT]), and Go in Eq. for various composite structures. O ur strategy in achieving 
this goal, which was strongly influenced by the w o rk of de Ge nnes and Okumura ( de Gennea . 1990l 2000; 



[5k umura and de Genn"el . l200ll:[OkumuraLl2002allbl . l2003Ll2005l) . is to write down coarse-grained linear elas 



tic energy functions and viscoelastic dissipation functions for biological composites of various structures, 
and to use available quasi-static crack solutions to estimate in a perturbative manner the small-velocity 
linear viscoelastic contribution to the fracture energy. We then use "matched asymptotics" considerations 
to show how this analysis can be extended to a wider range of crack velocities. We focus on the scaling 
properties of these quantities, i.e. we systematically neglect pre-factors of order unity, and consider three 
types of anisotropic structures: (i) Liquid-crystal-like composites (ii) Stratified composites (iii) Staggered 
composites, for different crack orientations with respect to the internal structure. 

In ordinary isotropic viscoelastic fracture (the details are provided i n sect ion 12. 2[) £ scales with a micro- 
scopic cutoff length (usually termed the "process zone size" (Broberg, 1999h ). E scales with the isotropic 



elastic modulus and £ scales with a macroscopic (geometric) cutoff length, the crack's length L for the frac- 
ture configuration considered here, and the functions g(-) and w(-) are of order unity. This is the hallmark of 
isotropic fracture: large scales elastic energy is dissipated at the small scales near the tip of a crack. Our re- 
sults, summarized in Tables 1-3, show that the presence of anisotropic nano-structures introduces additional 
lengthscales (determined, for example, by the ratio of the elastic moduli of the soft and hard constituent 
materials or by the aspect ratio of the plate-like mineral crystals) gives rise to different scaling behavior 



2 



as compared to isotropic fracture. In addition, we show that anisotropy may have some implications for 
fracture criteria. 

The structure of this paper is as follows. In Section [5] we describe the general procedure we adopt and 
apply it to isotropic viscoelastic fracture. In Section [3] we consider fracture in liquid-crystal-like structures. 
In Section 2] we consider fracture in stratified (layered) composite, for two crack orientations (parallel and 
perpendicular to the layers), while in Section [3] we consider fracture in staggered composites. Section [5] 
briefly discusses anisotropy effects and their relevance for fracture criteria. In Section [7J we show how to 
extend the perturbative approach to a wider range of finite crack velocities. Section [5] offers a summary and 
some concluding remarks. 



2. General procedure and application to isotropic viscoelastic fracture 

2.1. General procedure 

In order to calculate Go and G V i S for various nano-structures we describe first the general procedure we 
follow. In this Section we focus on a perturbative approach, while in Section [7] we show how to extend the 
analysis to the non-perturbative, finite v, regime. As str e ssed above, we have enormously benefitted from 
the pa pe rs of de Gen nes and Okumura (|de Germed . fl99fl l200ft lOkumura and de Gennesl . l200ll : lOkumural 



2002allbL 120031 120051) and in many ways the present contribution is a further development of their work 



The starting point of our procedure is a nano-mechanical model that incorporates the salient features of 
a given nano-structure into a continuum, coarse-grained, description of the effective viscoelastic response 
of a composite material. Specifically, this crucial step results in a coarse-grained elastic energy density 
function f(Vu,{Ei}) and a coarse-grained dissipation function i?(Vii, {rji}), where u is the coarse-grained 
displacement field and {Ei,rji\ is a set of elastic moduli and viscosity coefficients (respectively) associated 
with the different constituent materials in a given nano-composite. Note that the dimension of / is energy 
per unit volume and of R is energy per unit length per unit time. 

The next step in the perturbative approach is to derive an equation of motion for u by looking for the 
stationary variation of /(V«, {Ei}) with respect to u and neglecting i?(V«, {Vi})- Then one should solve 
the equation for u in the presence of a crack, i.e. a line that cannot support tensile and shear traction, under 
the external loading conditions. In many cases scaling arguments and matched asymptotics considerations 
can be useful in obtaining the scaling properties of u in different regions in space. 

The final step is to calculate Go and G V i S as follows. The two-dimensional crack configuration we 
described above is sketched in Fig. [1] The crack is assumed to propagate at a velocity v which is much 
smaller than the speed of sound. Recall that the presence of a crack is expressed by the usual mixed 
boundary conditions 

d±u± = on the crack , (4) 
uj_ = along the symmetry line, outside the crack . (5) 

The shaded area represents the typical strain distribution around the crack, where Aj_ is the scale in 
the direction perpendicular to the crack and Ay is the scale in the direction parallel to the crack. The 
relative magnitude of Ax and An will vary from problem to problem and in not necessarily as is shown 
schematically in the figure. The energy release rate Go is just the elastic strain energy density ahead of the 
crack tip multiplied by the spatial extent of the strain distribution in the direction perpendicular to crack 
propagation. Therefore, following Fig. [TJ we can immediately write it as 

Go ~ |jA ± , (6) 

where E± is the elastic modulus in the perpendicular direction. Comparing Eq. © with Eq. implies 
that E in the latter is E± and that Ax ~ Lg(£/L). The strain in the perpendicular direction, on a scale 
Ax, is easily estimated as 

A ± ~ E X ■ {l) 
3 




Figure 1: A schematic sketch of a crack of length L in an infinite medium (i.e. a medium of linear dimensions much larger 
than L) remotely loaded by a uniform tensile stress (Too. The shaded area represents the typical strain distribution around the 
crack, where Ax is the scale in the direction perpendicular to the crack and Ay is the scale in the direction parallel to the 
crack. The relative magnitude of Aj_ and Ay will vary from problem to problem and is not necessarily as shown schematically 
in the figure. 



Finally, the fracture viscous dissipation G V i S is given as 

G ms = . (8) 

v 

Under steady state propagation conditions we can replace time derivative in R(Vu, {r]i}) by space deriva- 
tives following dt — — where || denotes the direction of crack propagation. The procedure described here 
allows one to systematically calculate Go and G V i S in a large class of problems. It is important to note that 
this procedure is perturbative in nature since we calculate the static crack solution disregarding viscous de- 
formation, and use it later to estimate the dissipation arising from viscous deformation. A non-perturbative 
extension of this procedure is discussed in Section [7J 

2.2. Isotropic linear viscoelastic fracture 

In order to demonstrate how the procedure described above actually works and also to set a reference 
case for the anisotropic problems to be considered later, we start by discussing ordinary isotropic linear 
viscoelastic fracture. Isotropy implies that 

Ex,E\\~E, u±,u\\ ~ it, A_l,A||~L, (9) 

where u is the displacement in either the x or y directions and E is the isotropic elastic modulus. According 
to the first step in our procedure, we need to write down an expression for / and R. The elastic energy 
density has the form (scaling-wise) 

f~E(du)\ (10) 
where the spatial derivative d refers to either x or y. In addition, the viscous dissipation function reads 

R-r, J {dufrdrdO , (11) 

where rj is the viscosity, (r, 9) is a polar coordinates system whose origin is at the tip of the crack and the 
two-dimensional integration extends over the whole body. 

The next step is to derive an equation of motion for u by looking for the stationary variation of / with 
respect to u and to solve it in the presence of a crack. It is well known that the resulting equation is the 
Lame equation and that the solution at intermediate scales reads (|Brobersl Il999h 



which implies the well-known square-root singularity for the stress, a ~ (Too-J L/r. Intermediate scales mean 
scales between a microscopic inner cutoff lengthscale, say d c , where the linear theory fails (the so-called 
"process zone scale") and the outer scale L. As a consistency check we note that at r~L, we recover Eq. 
([7]) (recall that A±^L) and a^a^ as expected. 

The last step is to use Eqs. © and © to estimate Go and G V i S . Substituting Aj_ ~ L in Eq. ©, we 
immediately obtain 



Go ~ , (13) 



E 



We now substitute Eq. (fT2"j) in Eq. (0) to obtain 



S AWIW ~ E » x AW* ~ ^ . r /; * ~ Go « r (1 - i) ~ 



where <i c <§C L and t = j]/E. 

Comparing Eqs. (TT3|) and fj 14[) with Eqs. ((2J and ([3]) we immediately observe that w^O(l) and g~0(l), 
i.e. in isotropic fracture the viscous dissipation is controlled by a microscopic lengthscale l~rf c and quasi- 
static energy release rate is controlled by a macroscopic lengthscale t~L, This is the hallmark of isotropic 
fracture: energ y is released from large scales and is dissipated at the small scales near the tip of a crack 



(Brobcrg, 199 



3. Liquid-Crystal-like structures 

The structure of smectic liquid cr ystals is somewha t reminiscent of the structure of some biological 
composites, as was previously noted in lOkumural (|2002bh . We therefore start our discussion of anisotropic 
composite structures by considering liquid-crystal-like structures. By that we mean a structure that is 
composed of ordered layers of width d that lie along, say, the x-direction and can deform along the y-direction, 
by both stretching and bending. No elastic deformation along the x-direction takes place. Furthermore, 
gradients of the displacement rate (material velocity) give rise to a visc ous response. Such a material is 



characterized by the following elastic energy functional ( de Gennes . 1990h 



f ~ E{dyUy) 2 +Ed 2 (d XX Uyf 



(15) 



Here u y (x,t) is the displacement, E quantifies the elastic stiffness in the y-direction and Ed 2 is associated 
with the bending stiffness of the layers. The viscous dissipation function of such a material takes the form 



R 



[(a 



x J 



(dyiLy) 2 ] dxdy , 



(16) 



where 77 is the viscous response coefficie nt. Th i s mat erial response can approximately describe systems such 
as lamellar phases of block copolymers (jKatol . 120021) . but strictly not bone-like materials. Nevertheless, we 
believe that this example is very instructive and relevant in the presen t context. 

The stationary variation of / in Eq. fj 15[) with respect to u y reads ( de GennesI 1990) 







yyu,y 



d 2 d. 



xxxx U"y 



u qJ = 



(17) 



which is of course valid on scales larger than d. The crack is parallel to the layers, i.e. located along 
the cc-direction. The competition between stretching and bending (i.e. the different order of the spatial 
derivatives) in Eq. (|17[) immediately implies anisotropic scaling. Specifically, we observe that 



At - dA„ 



(18) 



where A x is Am of Fig. [TJand A y is Aj^ of Fig. [TJ Since we expect A X ^L ^> d, we immediately deduce that 



A„ 



L' 2 



(19) 



Substituting this result in Eq. ([3]) and identifying E±_ = E yields 



Go~^, (20) 



which is identical to the result previously derived in Ide Gennes (Il990h . We first observe that Go for liquid 



crystal-like structures contains the microscopic scale d. Furthermore, it is a factor L/d^$>l larger than the 
isotropic result of Eq. (fl"3|) . 

In order t o calculate G„i K w e need the solution of Eq. (|T7| in the presence of a crack. This problem was 



considered in de Gennes ( 1990t ). where it was found that 



u ^ x \Wv)- (21) 

Focusing on the outer scale of the problem, i.e. x^A x ^L and y^ A y ~L 2 /d, Eq. ([7]) immediately tells us 
that in fact 

a °° L xh (jL}\ . (22) 



Ed Vv^y 

To estimate G v i s , we first note that 

R.~i]v 2 J [(d xx u y ) 2 + {d xy Uy) 2 ] dxdy — rjv 2 J (d xx u y ) 2 dxdy , (23) 
because A y ^> A^, (i.e. gradients in the x direction are much larger than in the y direction). Then, we obtain 

G vis = — ~ i]v [ (d xx u y ) 2 dxdy ~ bt ^ ( [ —^^-dxdy , (24) 



Ed 2 J J dy 

where k(-) is related to h(-) and its derivatives, and r — rj/E. The function fc(-) in the integrand is finite in 
the limit y — > 0. Therefore, the integrand is characterized by an integrable singularity that scales as y~ x l 2 
(easily seen by introducing an auxiliary variable x = x/^/dy) and the integral is dominated by the upper 
limits of integration, which satisfy the scaling relations discussed above. Finally, performing the integration 
we obtain 

G v is ^ Go — — ■ (25) 
d d 

Comparing this result to Eq. @ we observe that d c ~ d, w ~ L/d and I ~ L. This suggests that G V i S 
is affected by the large scale L. It is remarkable that crack dissipation is controlled by the outer scale of 
a fracture problem. Wc believe that this unusual result derives from the infinite anisotropy of the liquid- 
crystal-like structures; those structures posses no elastic stiffness at all in the direction parallel to the layers. 
Finally, we note the if we interpret d c in Eq. (|14p as the width the layers d, which is rather artificial, we can 
say that G V i S for liquid-crystal- like structures is a factor L/(i> 1 larger than the result for isotropic linear 
viscoelastic fracture. 



4. Stratified composites 

Here we consider stratified (layered) biological nano-composites, see Fig. [2] These structures are periodic 
in the y-direction and are translationally invariant in the ir-direction. They are composed of hard layers 
of width d and elastic modulus Eh and soft layers of width d and elastic modulus E s <^Eh. The effective 
Hooke's law for this structure can be calculated exactly, as presented in detail in Appendix |Appcndix A| 
For a mineral volume fraction of the order of 1/2 the result reads 

/ - E h {d x u x ) 2 + E s {dyu y f + E s {d x u x ){d y u y ) + E s {d x u y + d y u x ) 2 + E h d 2 (d xx u y ) 2 , (26) 



G 



Figure 2: A schematic sketch of a stratified structure composed of repeated parallel layers of width d. The darker layers are 
elastically hard, with isotropic stiffness E^ and width d, and the brighter layers are elastically soft, with isotropic stiffness 
E a <C Eh and width d. 



where a bending term, identical to the one that ap peared for t he liqu id-crystal-like structures, was included. 
This expression is identical to the one proposed in lOkumural (l2f)02bh . 

In order to obtain the viscous dissipation function, we assume that the hard material is purely elastic 
and that the soft material is linear viscoelastic with a viscous response coefficient rj s . Therefore, we obtain 



R ~ Vs / [(dyiiy) 2 + (d x u x )(dyUy) + (d x u y + d y u x ) 2 ] dxdy 



(27) 



Unlike the liquid-crystal-like structures, stratified composite structures have elastic stiffness in both the 
x and the y directions. Therefore, we should (at least) distinguish between cracks that are parallel and 
perpendicular to the layers. 

4-1. Parallel cracks 

We consider a crack located along the x-direction, i.e. parallel to the layers. The presence of a new small 
parameter 

F 

(28) 



implies a new lengthscale in the problem 



We focus on the limit 



A = -= > d 



d < A < L . 



(29) 



(30) 



As was shown in lOkumural (|2002bf ) , in this case the elastic functional of Eq. (f2l)|) cannot be approximated 
by a single expression for all relevant scales. We theref ore consider separ ately large scales, x^X, and small 
scales, :£<§; A. Consider first the large scales. Following Okumura ( 2002bl ). the elastic functional of Eq. (|26[) 
can be approximated by 



f ~ E s (d x u y y + E s (d y u y y 



(31) 



This implies isotropic, mode Ill-like, fracture in a soft material. We therefore use the results of isotropic 
fracture, now with a lower cutoff scale A 



G 



Go 

(0 



4. 



Go- 



(32) 
(33) 



where the superscript (I) denotes large scales. 

For the small scales, the elastic functional of Eq. (J^HJ) can be approximated by 

/ ~ E s {d y u y f + E h d 2 {d xx u y f . (34) 

To proceed, we employ "matched asymptotics" considerations, i.e. we determine the properties of the 
"inner" (small scales) solution by demanding that it smoothly crosses over to the "outer" (larger scales) 
solution at a scale x~ A. Therefore, the inner problem is characterized by a loading cr£^ ~a ca ^fL~/~\, with 

E®=E„, (Af)) 2 ~ AAjf), A0~A«=iA. (35) 

Here the superscript (s) denotes small scales, i.e. all of the quantities refer to the region near the tip where 
the energy functional in Eq. (|3"4")l dominates the elastic behavior and not to the shaded region in the global 
problem sketched in Fig. [1] 

Go for the small scales must be identical to that of the large scales as no dissipation takes place in the 
intermediate region. Using Eq. (J7J), we obtain 

^ ^/i(^=] . ism 



y E s V A VVAy 
where we used x~A. x ~\. We then note that 

R~ Vs J [(dyUy f + {d x u x ){d y u y ) + (d x u y + d y u x ) 2 ] dxdy ~ E s tv 2 J ' (d xy u y ) 2 dxdy , (37) 

because A y -C A. x (i.e. gradients in the y direction are much larger than in the x direction). For G^ s we 
obtain 



Git ~ E srv I (d xy u y ) 2 dxdy ~ E s rv^- I _ I "^ dyJ dxdy , (38) 



e 2 ^ Jx^vxdJy^d y 2 

where fc(-) is related to h(-) and its derivatives. The function fc(-) in the integrand is finite in the limit 
y — > 0. Therefore, the integrand is characterized by a non-integrable singularity that scales as y~ 3 / 2 (easily 
seen by introducing an auxiliary variable x = x/^/dy) and the integral is dominated by the lower limits 
of integration, which satisfy the scaling relations discussed above. Finally, performing the integration we 
obtain 



G V is ~ Go— 7= • (39) 
V« A 



The overall G V i S is then given as 



G vis - G„} 3 + G v f s — G — ( 1 + y ^ J — ^°^7^ ' ( 40 ) 
It is interesting to note that the result in Eq. (|4"U|) can be written as 

Gyis ~ G ^ , (41) 
an 

where d\\ is the smallest cutoff length in the crack propagation direction. In the present case the crack 

propagates in the x-direction and we have d\\ =d x ^\fd\. As we will see below, this result is rather generic 
as long as the dissipation integral is dominated by the lower limits of integration. 



8 



4-2. Perpendicular cracks 

Here we consider a c r ack lo cated along the y-direction, i.e. perpendicular to the layers. Following 



Okumura and de Genned (|200lh . in this case the elastic functional of Eq. (f^5|) can be approximated as 

f ~ E h {d x u x ) 2 + E s (d y u x ) 2 . (42) 
Minimizing the elastic energy with respect to u x , we obtain the following equation 

d XX U X + 8yyU x = , (43) 

wher e y = t~ x l 2 y (recall that e = E s / Eh <C 1). This is analogous to a mode-Ill crack problem ( Brober d , 
19991) . 



Note that for the crack orientation considered here Aj_ =A X , A|| = A y and E± = Ey l . Since we expect 
AyC^L, we immediately obtain 

A x ~ e-^L » A y . (44) 

Substituting in Eq. ([5]), we obtain 

This result shows that the energy release rate in this case is significantly enhanced as compared to fracture 
in an isotropic material with stiffness Eh- 

We consider now Eq. ([7]) which tells us that the crack opening on a scale y ~ L reads 

Ux „ £^ e -V2 £ . (46) 

Furthermore, from Eq. (|43[) we know that the opening displacement should satisfy ordinary fracture scaling 
in the coordinates (x, y), i.e. 

u x ~(x 2 +f)^ = (x 2 + e-V) 1/4 ■ (47) 

Therefore, we conclude that 



Ux-^e-^vT^+e-V) 1 /*- (48) 
We can now use the above expression for to calculate the viscous dissipation. We first note that 

R ~ Vs J [{dyUy) 2 + {d x u x )(d y Uy) + (d x u y + d y u x ) 2 ] dxdy ~ E s tv 2 J (d yy u x ) 2 dxdy , (49) 

because A y <^ A x (i.e. gradients in the y direction are much larger than in the x direction) and dt = —vd y . 
Therefore, we obtain 



G« s = — ~ ^ s Tf / (d yy u x ) dxdy ~ E s ra °° . / / -p—: — dxdy , (50) 



F 2 ^l/2 / - ,„ / fl/2,/3 



where we used the lower limits of integration alone because of the small scales divergence. The integral 
alone yields (ed)^ 1 , which implies that 

G ^-E s rv^-^-^G -. (51) 

Note that this result agrees with Eq. (|4"Tj) . where this time we have dn —d y ^d. 
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Figure 3: A schematic sketch of a staggered array of hard platelets (darker) embedded within a soft matrix (brighter). The 
width of the platelets is d and their length is pd, where p is a dimensionless aspect ratio. The length of the horizontal (along 
x) gap between platelets in the same layer, filled with a soft material, is a. The width of the infinitely long soft layers is also d. 



5. Staggered composites 

Micro-structural studies ( Fratzl , 2007 ; Dunlop and Fratzl , 2010t ) have demonstrated that bone-like ma- 



terials feature a staggered arrangement of hard mineral platelets at the nano-scale, see Fig. [3] The width 
of the platelets is d and their length is pd, where p is a dimensionless aspect ratio. The length of the 
horizontal (along x) gap between platelets in the same layer, filled with a soft material, is a. The width of 
the infinitely long soft layers is also d. Our goal here is to find out whether the staggered structure gives 
rise to a different viscous fracture dissipation as compared to stratified structures. The structural difference 
between the staggered array and the stratified one is that in the former the hard mineral platelets are of 
finite length, while in the latter they are infinite. This implies that the resulting stiffness in cc-direction may 
be different. Therefore, denoting the elastic stiffness in this direction by E%/\ we can write 

/ ~ E e J f {d x u x f + E s (d y u y ) 2 + E s {d x u x ){dyu y ) + E s (d x u y + d y u x f + E h S '(fl ^ xx u y f . (52) 



In order to est i mate E%** , we ad opt the Tension-Shear-Chain (TSC) nano-mechanical model (jArzt et al 



2003; Ji and Gao, 2004; Gao, 2006). The major assumption of this model is that tensile stresses are trans- 



ferred between the hard mineral crystals mainly by shear stresses in the soft matrix. In other words, the 
tensile stress in the gap of length a is assumed to be negligible. This assumption is valid in the limit of large 
stiffness contrast £<§;!, large aspect ratio of the hard platele ts, p^$>l, and a^pd. In this limit, the model 



predicts the following expression for E x '* (jArzt et al 



l piateiei 
fl|2Q03j) 



+ , 1 (53) 



E e J f ~ E s <t> 2 p 2 <t>E h ' 

where v s is the soft material Poisson's ratio and <j> is the volume fraction of the hard material. We focus 
on 0~ 1/2. The first term on the right hand side of Eq. (j53l) represents the shear deformation in the soft 
material between the hard platelets, while the second one the tensile deformation inside the platelets. 

We are mainly interested in the regime where E s <C E s p 2 < Eh, i.e. when the shear deformation of the 
soft material is important. Under this condition we can write 

E e Jf ~ E s p 2 . (54) 

The corresponding dissipation function reads 

R~ r/ s J [p 2 {d x u x ) 2 + {dyiiy) 2 + {d x u x ){d y u y ) + {d x u y + d y u x ) 2 ] dxdy . (55) 

Comparing this expression to Eq. (1271) for stratified materials, we immediately identify the additional term 
proportional to p 2 (d x u x ) 2 . The physics behind this new contribution is clear; it represents the additional 
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viscous dissipation due to shear deformation of the soft matrix during tensile deformation in the x-direction. 
Since p3>l is a large number, one may expect a qualitative change in the fracture dissipation. 

In order to check if this is indeed the case, we first note that since for cracks parallel to the platelets (i.e. 
along the x-direction) the deformation is dominated by u y , the results are identical to those for stratified 
structures. For cracks perpendicular to the platelets (i.e. in the y-direction), we compare Eqs. (|52p and 
(I5~il) with Eq. (|2"B1) . which immediately suggests that we should simply identify Eh is the latter with E s p 2 , 
and use all of the results obtained before for perpendicular cracks in stratified structures. In particular, we 
identify e of Eq. (12"51) with p~ 2 and obtain 

E ± =E s p 2 , A x ~pA y , A y ~L, A x ~ pL 3> A v , (56) 

together with 

Go~^, u x ~-^VZ(* 2 +pV) 1/4 . (57) 
E s p E s p i l I 



We can now use the above expression for u x to calculate the viscous dissipation. We use Eq. ((55)) to 
obtain 



R ~ r] s J [p 2 {d x u x ) 2 + (dyiiy) 2 + (d x u x )(d y Uy) + (d x u y + d y ii x ) 2 ] dxdy 

~ E s tv 2 / [p 2 (d xy u x ) 2 + (d yy u x ) 2 ] dxdy , (58) 



The second term in the above expression must result in the same contribution as in Eqs. (1501) and (I51[) . 
We therefore focus on the first term, which is new and, as mentioned above, may possibly lead to enhanced 
dissipation. This term yields 

E s rvp 2 [{d xy u x ) 2 dxdy~E s Tvp 2 ^f^ f f k ^ PV ) dxdy . (59) 

J ^sP Jx~pdJy~d PV 

The integral alone yields d~ l . Therefore, this contribution to the viscous dissipation reads 

E s rvp 2 J (d xy u x ) 2 dxdy ~ E s Tv^- d ~ G ^ , (60) 

which is identical to the contribution from the second term. We therefore conclude, somewhat surprisingly, 
that the viscous dissipation accompanied perpendicular crack propagation in staggered structures are en- 
hanced compared to stratified structures only by a multiplicative factor of order unity, i.e. scaling-wise it 
remains unchanged 

TV 

G V is ^ Go— . (61) 

The origin of this unexpected result is that while the dissipation associated with the deformation in the 
x-direction in Eq. (|59[) is indeed enhanced by a large factor p 2 , the displacement u x in Eq. (|57[) is reduced 
in a way that precisely cancels out the p dependence of the viscous dissipation. 



6. Anisotropy and fracture criteria 

Our goal here is to demonstrate, through a brief example, that the anisotropy of the structure and 
the constitutive response may have implications for fracture criteria. The example we focus on is a crack 
perpendicular to the layers in a stratified structure (i.e. located along the y-direction). Using standard 
fracture mechanics, Eq. (j43j) immediately implies 

e~ 1/4 goo%/£ (f . \ 
Uxx ~ (x 2 +e~V) 1/4 ' 
11 



The prefactor e 1 / 4 ensures that a xx approaches tJoo on a scale y ~ L, as expected. Focusing on the crack 
symmetry line, x = 0, we obtain 

cr xx (x = 0,y) ~ . (63) 

This result shows that the near tip st ress is unchanged as compared to isotropic fracture (i.e. the effective 



stress intensity factor ( Broberg , 19991 ) remains K e ^ ~ a ao \fV)- However, the energy release rate for this 



configuration, already calculated in Eq. (|4"5j) , is significantly enhanced as compared to isotropic fracture (by 
a factor e -1 / 2 ). This shows that the standard isotropic relation between the stress intensity factor and the 
energy release rate may not be valid in anisotropic situations. In addition, this comment suggests that some 
caution should be taken in relating a fracture toughness criterion (i.e. a critical stress intensity factor K c ) 
and a fracture energy criterion (i.e. a critical energy release rate G c ) in anisotropic materials. 

To further strengthen the point about such anisotropic effects, we briefly consider a related fracture 
problem in which a long crack (again perpendicular to the layers in a stratified structure) is located on the 
symmetry line of an infinite strip of width W and loaded by a fixed-grip boundary condition at the lateral 
edges of the strip. In this case the relevant macroscopic lengthscale in the proble m is W and not the crack 
length (which is assumed to be much larger). This problem was considered in Ok umura and de Gennesl 
(2001). Global energy balance during crack propagation immediately implies that 



a 2 W 

Go ~ -|- , (64) 

where gpo here is the homogeneous str ess far ahead of the crack tip (even though there is no real "infinity" 
here). [(5k umura and de Genned (|200ll) found that 



..(..O,,)-^? and ^M- r y , (65) 

which is indeed consistent with Eq. (|64[l through Go ~ cr xx u x . Comparing Eqs. (|64p and (|65[) with Eqs. 
(|45p and (|63p demonstrates that fracture has different properties in strip geometry and in infinite medium 
for anisotropic materials. Therefore, again, some caution should be taken in interpreting fracture criteria. 



7. The higher velocities regime 

The results presented up to now are restricted to small crack velocities, i.e. they were derived in the 
framework of a perturbation theory. For example, in the case of stratified and staggered composites, the 
relevant velocities range is vr/dti < 1 (where dn may vary from one problem to another). It is important to 
note that strictly speaking the perturbative approach is valid for vt /d\\ <C 1, which is the small parameter 
in the perturbative expansion, but as usual the results provide sensible estimates up to vr/d\\ ~C(1). Our 
goal here is to demonstrate, through an example, how one can obtain G V i S (v) at finite velocities without 
solving the equations of motion and evaluating the dissipation integral. 

To elucidate the procedure we adopt, we first discuss isotropic viscoelastic fracture. We start with the 
perturbative result, i.e. Eq. (|14[) . which is presented here again for completeness 

G ms ~ G ^ . (66) 

As explained above, it is valid for vt / d c < 1. Our aim is to extend this result to higher velocities, vr/d c > 1. 
We first note that at higher velocities rate-dependent effects may emerge. That is, regions at different 
distances from the crack tip, r>d c , may experience qualitatively different strain rates e through the relation 

e = v/r . (67) 

In particular, with increasing velocity the region near the tip of a crack may experience strain rates £3>t _1 
which may result in a modified material response. In fact, we expect disordered (e.g. polymeric) materials 
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to exhibit significant stiffe ning when deform ed at rates much larger than their typical mechanical relaxation 
rate. Therefore, followinglde Genned(|l996l) . we assume that the material under consideration is viscoelastic, 
characterized by E s and r] s , at small strain rates and elastic, characterized by Eh^>E a , at large strain rates. 
For sufficiently large crack velocities a region characte rized by Eh devel ops around the crack tip. Under 
these conditions G V i S has been shown to take the form ( de Gennes . 19961) 



^vis — "0 7T 



(68) 



The question is then how the perturbative, small v, result of Eq. (|66| is smoothly connected to the large v 
re gime of Eq. (|6"8|). Th i s prob lem was addressed and systematically solved, in a broader viscoelastic context, 



Persson and Brenerl (12005 ). We fo llow here the approach of the latter paper. The important observation 



made in (jPersson and Brenerl 120051 ) is that the cutoff scale d c is in fact a dynamic quantity that grows with 
v, d c (v), and hence we should interpret all previous appearances of d c as d c {v — 0). Furthermore, it was 
assumed that d c (v) evolves such that the stress level at that distance from the tip is constant, i.e. that there 
ex ists an ^-independen t yield stress or some analogous stress quantity. With this assumption, it was shown 



tyii 

Persson and Brenerl (|2005[ ) that 



Gtotiv) d c (v) 



4(0) 



(69) 



where the dissipative contribution emerging from the region r < d c {v) (plastic and process zone) was neglected 
as compared to the viscoelastic contribution emerging from the region r > d c (v). While we do not think 
these assumptions are universally valid, we do believe they provide a possible framework to make progress 
and gain insight into the fracture dissipation in this class of problems. 

With Eqs. (|66p . ((68|) and (j69|) at hand we can now complete the calculation using "matched asymptotics" 
considerations. We first require that Eq. (|66[) is approached in the limit v — > and hence replace d c there 
by d c {v). Eliminating d c {v) between the resulting equation and Eq. (|69[) . using Gtot — G V i S and solving for 
G V i S , we obtain 




G r 



Since this result should smoothly connect between Eqs. 
range of velocities 



(70) 



(|6"6")l and (|68p . it must be valid in the following 



— < v < 

T 



E h 



dc 

T 



(71) 



This provides a description of G V j S (v) bey ond perturbation theory and over a wide range of crack velocities, 
in accord with iPersson and Brenerl (|2005l ). We stress that inertia, as in the rest of this paper, is negligible 
is the present context. 

We now demonstrate how these ideas are applied to anisotropic viscoelastic fracture of biological com- 
posites. We focus on large cracks that propagate parallel to the layers in a stratified structure and assume 
that the soft constituent material has a response similar to the material considered above in the isotropic 
case. The hard material, as elsewhere in this paper, is purely elastic with a modulus Eh at all strain rates. 
The perturbative result of Eqs. P0")) and pTj) reads 



Gq- 



Gq- 



<d\ 



(72) 



In this case the small scale e?|| =d x is related to the characteristic small scale in the y-direction, d y , by the 
relation d\^ Xd y . For small velocities, v< ^/\d/r= ^\/d(d/r) (i.e. in the perturbative regime), we know 



that d x ^vXd and d y ~d. For larger velocities we assume, as in the isotropic case, that Eq. (f72"|) still holds 
with d\\ =d x (v) defined such that the stress at this scale remains w-independent. This leads to 



G vis (v) d y (v) d 2 x (v) 



Go 



dy(0) 
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(73) 



Using the same procedure as before (i.e. eliminating d x [v) etc.) we obtain 



\ 2/3 

G ms ~ Go ( ) . (74) 



XdJ 

This scaling holds until d x ~d y ~\, which corresponds to vr~A 2 /d, where fracture becomes isotropic. The 
latter is easily obtained by solving for d x (v) using Eqs. ((73|) and (|74| . Therefore, Eq. (|74[) is valid in the 
following range of velocities 

[Xd X 2 d 

VdT <V< ^T- (75) 

For vr>\ 2 /d fracture is isotropic and is controlled by the soft constituent material, cf. Eq. (13"TT) . Therefore 
the above results for isotropic fracture hold and we should determine the range of validity of each expression 
using again "matched asymptotics" considerations. In particular, for vt > X 2 /d we obtain the isotropic 
result, cf. Eq. JTOJ, 

/ VT 

G v is ~ ~d ' 

which is smoothly connected to Eq. (174")) at ur~ X 2 /d. This result is valid until vr/d~ (Eh/E s ) 2 , where it 
crosses over to (cf. Eq. ([55])) 

G„t S = Go^- . (77) 

The complete results for the example worked out in detail here are summarized in Table 3. It is easily 
confirmed, by direct substitution, that the result in each regime smoothly crosses over to the result at the 
next regime at yet higher velocities. Therefore, we have demonstrated how one can obtain the behavior of 
G V i S (v) for an anisotropic structure for a wide range of velocities, going beyond the theory of perturbation. 



8. Summary and conclusions 

In this paper we calculated the scaling properties of the energy release rate and the fracture viscous 
dissipation for various anisotropic composite structures. Our results are summarized in Tables 1-3. The 
dominant lengthscales for energy release and viscoelastic dissipation can be identified together with the 
non-trivial scaling behavior emerging from the various anisotropic composite structures. 



Table 1: Go and G v i a /Go for slow cracks (perturbative regime) parallel to the hard constituent material in the various structures 
considered in this paper. Recall that \ = dyfE^JE s . 



Slow cracks parallel to the hard constituent material 




Liquid-crystal- like 


Stratified composites 


Staggered composites 


Go 


Ed 


E, 


E. 


G v is /Go 


VT L 

d d 


VT 

Vd\ 


VT 

VdX 



Table 2: Go and G v i a /Go for slow cracks (perturbative regime) perpendicular to the hard constituent material in the various 
structures considered in this paper. Recall that p is the aspect ratio of the hard platelets in the staggered structure. 



Slow cracks perpendicular to the hard constituent material 




Liquid-crystal- like 


Stratified composites 


Staggered composites 


Go 




s/E,E h 


E 3 p 


Gvis /Gq 




VT 

A 


VT 
d 
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Table 3: G v i s /Go for cracks parallel to the hard layers in stratified composites for a wide range of velocities, i.e. in the 
non-perturbative regime. Recall that \ = d\J E^/ E s 



Cracks parallel to the hard layers in stratified composites 


Velocity range 


V d r 


\ d , A 2 d 

\d- <v< d^7 


>?r < < (E 1 X l d 
W d^ u ^\E 3 J t 


->{i) Z± T 


G V i S /Go 


VT 

-JYd 


( VT \ 


PUt~ 
V d 


E s 



We would like to highlight the qualitative difference between the viscoelastic dissipation in liquid-crystal- 
like composites and the other composites considered in this paper. The dissipation integral for the liquid- 
crystal- like composites, Eq. (|24l) . is dominated by the upper limits of integration, which is a result of the 
integrable singularity of the integrand. The consequence is that viscoelastic dissipation in this case comes 
from a region whose size is controlled by a macroscopic, extrinsic, material- independent, lengthscale. In 
contrast, the viscoelastic dissipation of other composites discussed here, even when involves different scales 
as compared to isotropic fracture, is always controlled by intrinsic, material-dependent, lengthscales. This 
derives from the appearance of a non-integrable singularity in the dissipation integral, which is a generic 
feature of ordinary fracture. We believe that this qualitative difference emerges from the fact that the 
anisotropy of the liquid-crystal-like composites is infinite, i.e. the direction parallel to the layers exhibits 
no elastic response at all (i.e. it features a liquid response in this direction), while in all other cases the 
anisotropy may be (very) large, but always finite. 

We believe that our results may have both theoretical and practical merit. From a theoretical point 
of view, our results show that in the presence of heterogeneous structures and variability in the local 
mechanical properties, the linear elastic energy release rate and the fracture viscous dissipation may exhibit 
different properties as compared to their isotropic fracture counterparts. In addition, our work further 
demonstrates the power and potential of continuum approaches in understanding biological (and other) 
composite materials. Note, however, that we cannot exclude other important effects emerging from scales 
that are not resolved in the coarse-grained approach. 

From a more practical point of view, our work highlights the possible importance of viscous dissipation as 
a toughening mechanism in biological composites. This may be at least partially supported by the difference 
between the fracture toughness of hydrated and dehydrated biological composites. We suspect that this 
velocity-dependent fracture toughness enhanceme nt may be quantitatively important, in addition to the 
fracture toughness associated with crack extension ( Nalla et al. . 2004 : Launev et al. . 2010). However, careful 



experiments are needed in order to test our suggestion. Such experiments require controlling the applied Gtot 
and accurately measuring the critical conditions at crack initiation Gq = G c , in order to calculate G V is/Go- 
Such measurements also involve monitoring the crack velocity v and independently extracting material 
parameters such as t and E s / E^. At present, we could not find extensive and systematic experimental data 
on such quantities. 

Finally, our results demonstrate the relevance of anisotropy to the interpretation of different fracture 
criteria in biological composites. This may be useful for interpreting and characterizing the fracture-related 
material properties of these composites. In this context, it is important to note that we foc ussed on strong 



aniso tropy (either in elastic or fracture propert ies) , which is indeed observed at small s c ales rtPeterlik et al 



2006), but usually not on macroscopic ones, cf. iBehiri and Bonfieldl (|l989h : iNalla et all (|2003l ). The latter is 



understood as resultin g from addit i onal, larger scale, s t ructu res that are common in multi-level hierarchial 
biological composites ( Fratzl , 2007 ; Dunlop and Fratzl . 2010l ). This may be the subject of a future investi- 



gation. 
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Appendix A. Exact coarse-grained elastic energy for stratified structures 

Here we derive the exact coarse-grained linear elastic response of stratified (layered) structures. The 
hard material has a Young's modulus Eh and a Poisson ratio i//,, and a thickness dh- The soft material 
is characterized by E s and v s , and a thickness d s . The volume fraction of the hard material is given by 
</> = dh/(dh + d s ). Let us first derive the effective Hooke's law for this structure. We treat the hard and 
soft materials as isotropic and linear elastic, and the interfaces between them as perfect. We denote the 
periodicity direct ion by y and the pe r pendi cular direction by x, and assume plane stress conditions. We 
therefore obtain ( Landau and Lifshitz . 19861) 

J h ) _ „ Jh) Js) _ (a) Jh) _ „ Jh) (a) _ (a) 

Jh) _ U VV "hUxx J s ) _ Uyy v s" xx Jh) _ u h"yy J s ) _ "xx ^s"yy , . ^ 

yy ~ E h ' yy ~ e s ' xx ~ E h ' xx ~ e s ' [ ' 



4h) Jh) = Oxy{l + Vh) = CTy h x \l + V h ) J s ) = ( S ) = CTxyjl + Vs) _ (7$ (I + 1S S ) 

--xy yx ^ ^ > xy yx ^ ^ 



Here we used the definition of the linear elastic strain tensor 



Jh,a) 



IfaM+djuM) (A.3) 



and the conservation of angular momentum in each material u x h y S ^ = u y l x ' s \ The superscripts (s) and (h) 
correspond here to the soft and hard layers, respectively. 

The boundary conditions demand that the displacement along the interface is the same in both materials 
(i.e. a perfect interface) and that the force across the interface is continuous. Since the interface is placed 
along the x-axis (constant y), these conditions translate into 

<J x u- x —u x u x , u x a y —u x tiy , u yx —u yx , °yy "yy • K 1 ^-^) 

Using angular momentum conservation we also obtain a X y = a x h y . Let us first consider the extensional 
components. The average o xx is given by cr xx = (l — <f)cr x b x +4>a xx . 
The condition e xx — e x s x — e x h J implies 



(h) _ (h) (a) _ (s) 

O xx "h"yy &xx ^s^yy 

Eh E s 

Solving the last two equations we obtain 

Oxx , / Eh 



(A.5) 



V s - Vh CTyy , 1 ,s (s) 

(s) _ f_ \E S / (h) _ ^ra _ I 1 — <P)Vxx 

T xx - E h \~4> ' ° XX 



(A.6) 



E s 

Substituting these expressions in 



(ft) , /-, _ ,N (a) _ " VhVxx) (1 - <P){(Jyy - VsOxx) _ (J xx ~ V s <Jyy 

-yy + U Phyy ^ + ^ , ?xx - ^ 



results in the extensional part of the effective Hooke's law, that is, a linear relation between e xx ,e yy and 
OxxT&yy For the shear part we have d x u y = d x u y s ^ = d x Uy . Therefore, 



E s [d x u y + d y u x s J E h (d x u y + dyui^ 
2{l + v s ) ~ 2(1 + u h ) 



• X y ■ 



(A., 
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Expressing d y u x ' in terms of a xy and using d y u x = (1 — <f>)d y Ux +4>d y u x , we obtain 

a xy . (A.9) 

Using these expressions, the elastic energy density EijUij/2 (excluding the bending energy) can be im- 
mediately calculated. Assuming Eh 3> E s and </> ~ 1/2, adding the bending energy and omitting prefactors 
of order unity, we obtain Eq. (|26[) in the text. 
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